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■ Abstract. The present lecture notes address three columns on which the Lagrangian 
perturbation approach to cosmological dynamics is based: 1. the formulation of a 
Lagrangian theory of self-gravitating flows in which the dynamics is described in terms 

^ of a single field variable; 2. the procedure, how to obtain the dynamics of Eulerian 

ly-j ' fields from the Lagrangian picture, and 3. a precise definition of a Newtonian cosmology 

. framework in which Lagrangian perturbation solutions can be studied. While the first is 

I a discussion of the basic equations obtained by transforming the Eulerian evolution and 

' field equations to the Lagrangian picture, the second exemplifies how the Lagrangian 

■ theory determines the evolution of Eulerian fields including kinematical variables like 
expansion, vorticity, as well as the shear and tidal tensors. The third column is based on 
a specification of initial and boundary conditions, and in particular on the identification 
of the average flow of an inhomogeneous cosmology with a "Hubble-flow" . Here, 

2 ' we also look at the limits of the Lagrangian perturbation approach as inferred from 

, comparisons with N-body simulations and illustrate some striking properties of the 

solutions. 
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1. Lagrangian Theory of Self— gravitating Flows 

The description of fluid motions in cosmology has been largely studied in an Eulerian 
coordinate system rectangular non-rotating frame in Euclidean space. Quite 

recently, it has become popular to study fluid motions in a Lagrangian coordinate system 
X, i.e., a curvilinear, possibly rotating frame in Euclidean space which is deflned such 
as to move with the fluid. Since the Lagrangian description has a number of advantages 
over the Eulerian one, and since this description enjoys many applications in the recent 
cosmology literature, it is important to elucidate in proper language the Lagrangian 
formalism. Since the lectures by Frangois Bouchet and Peter Coles (this volume) explore 
the fleld of recent applications of the Lagrangian perturbation theory, I here concentrate 
on the basic architecture of a Lagrangian theory of structure formation. I do this in 
Newtonian cosmology, the lecture by Sabino Matarrese (this volume) gives an extension 
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to the framework of General Relativity. Accordingly, I kept my reference list short, since 
more references may be found in the other lectures. 

That the Lagrangian approach is experiencing a revival in cosmology is good news; I 
consider it the natural frame to describe fluid motions, since this description is formally 
close to the mechanics of point particles. If you consult old textbooks on hydrodynamics, 
you will find that the Lagrangian picture was considered too complicated for practical 
purposes beyond problems with high symmetry, and therefore has not been pursued 
further. I hope that, after this and the related lectures, you will be convinced of the 
opposite. 

Let us start with the basic system of equations in Newtonian cosmology describing 
the motion of a pressureless fiuid in the gravitational field which is generated by its 
own density. We think at applications for a (dominating) coUisionless component in the 
Universe; the gravitational dynamics we describe is thought to act as an attractor for 
the baryonic matter component which is "lighted up" by physics not described by these 
equations. 

With this assumption the fiuid motion in Eulerian space is completely characterized 
by its velocity field v{x, t) and its density field g{x, t) > 0. The fluid has to obey the 
familiar evolution equations for these fields. 



where the gravitational field g{x, t) is constrained by the (Newtonian) field equations 



A denotes the cosmological constant. (Strictly speaking, g in eq. (la) is a force per unit 
inertial mass, whereas in eqs. (lc,d) g is the field strength associated with gravitational 
mass. That we set both equal is the content of Einstein's equivalence principle of inertial 
and gravitational mass.) 

Alternatively, we can write the eqs. (Ic) and (Id) in terms of a single Poisson 
equation for the gravitational potential, which we do not need in the following. 
Hereafter, we call the system (1) the Euler-Newton system. 

One important issue to learn about the Lagrangian treatment is the fact that both 
evolution equations (la) and (lb) can be integrated exactly in Lagrangian space, velocity 
and density will therefore not appear as dynamical variables later. To see this, we first 
look at the basic Lagrangian field variable which is the trajectory field of fluid elements, 
or the deformation field of the medium, respectively (Fig. la): 



dfV = -{v ■V)v + g , 
dtQ = -V • (gv) , 



(la) 
(lb) 



Vx^=0, 

V • ^ = A - AttGp ; 



(Ic) 
(Id) 



f =/(X,i) ; X:^f{X,to) , 



(2a) 
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where X denote the Lagrangian coordinates which label fluid elements, x are the 

— * 

positions of these elements in Eulerian space at the time t, and / is the trajectory 

— * 

of fluid elements for constant X. (Notice: the Eulerian positions x are here viewed 
not as independent variables (i.e., coordinates), but as dependent fields of Lagrangian 
coordinates; therefore, we employ the letter / for the sake of clarity. Independent 
variables are now {X,t) instead of {x,t).) 




Figure 1. The Newtonian spacetime projected onto an Eulerian plane is scetched. 
A fluid element sitting initially {t = to) at X moves along the trajectory / to the 
Eulerian position x at time t: this position vector is expressed in terms of Lagrangian 
(i.e. initial) coordinates which are constant along the trajectory: X = (Fig. la). 
The Jacobian J measures the volume deformation of fluid elements located a,t X; J 
drops from 1 to zero as the volume element degenerates into a surface element (piece 
of pancake), a line element (piece of filament), or a point (cluster). 



The velocity field v{x, t) of the fluid is the tangential field to this family of curves, 

v='f{X,t) , {2b) 
and the dot denotes the total (Lagrangian or convective) time-derivative along v: 

{...)' ^j^--^d, + v-V . (2c) 
(The dot commutes with Lagrangian differentiation.) 

In view of (2c) we recognize that the Eulerian evolution equation (la) can be written as 
V = g, which is the familiar definition of acceleration (here, the acceleration of a fluid 
element along its trajectory). Consequently, in view of (2b), we know the acceleration 
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field in Lagrangian space for any given trajectory: 



{2d) 



In otlier words, taking (2b) and (2d) as definitions of velocity and acceleration, 
respectively (like in point mechanics!), we see that equation (la) is automatically fulfilled 
as we switch to the Lagrangian description. 

A similar logic applies to the continuity equation (lb): the deformation of the 
medium is described by the Lagrangian deformation tensor (f^k), i-e., the tensor of 
first derivatives of the trajectory field with respect to Lagrangian coordinates f which 
measures how much the Eulerian positions deviate from their original (Lagrangian) 
positions. The volume of the deformed fluid element is measured by the determinant of 
this tensor, J{X, t) := det(/i,fc), and therefore its density must be inversely proportional 



of the Lagrangian coordinates in (2a) which coincide with the Eulerian ones at t = to-) 
You may verify (2e) by differentiation and by using the identity 



Again we conclude that, for any given f{X,t), wc obtain the exact expression for the 
density at the fluid element, i.e., in Lagrangian space. 

Let us now move to the question which Lagrangian evolution equations the field of 
trajectories has to obey. We obtain them by expressing the field equations (Ic) and 
(Id) (which are 4 hnear "constraint equations" for the acceleration field) in terms of 
Lagrangian coordinates which, as we shall see, yields 4 non-linear partial differential 
equations for g{X,t). For this purpose we need the inverse of the transformation (2a) 
from Lagrangian to Eulerian coordinates. 



While the deformation tensor is the Jacobian matrix of the transformation from Eulerian 
to Lagrangian coordinates, fi.k — Jik, the inverse of fi.k is the inverse Jacobian matrix 




(2/) 



X^h{x,t) ] h = f 



t 




Jab — 0'd{Jab)J ^ — -Z-:^ajk^bimfi\jfm\k ■ 



{2h) 



t Lagrangian differentiation is indicated throughout this lecture with a vertical slash to make a 
difference to differentiation with respect to Eulerian coordinates, denoted by commata. 
I Hereafter we adopt the summation convention. 



5 



Since the transformation part "looks" inconvenient during a first reading of papers on 
Lagrangian models, I here try to be as elementary as possible. 

We look at the two-dimensional case first, i.e., there are only two non-vanishing 
field components, e.g., gi and g2 f- We then get for the Eulerian derivatives of g: 

9l,l = 9l^lhl,l + 9l^2h2,l , fl'1,2 = fl'l,l/il,2 + 5'l|2/i2,2 , 

52,1 = fi'2|l/il,l + 5'2|2/i2,l , 5'2,2 = 5'2|l/il,2 + 5'2|2/i2,2 • (3a) 

Inserting the inverse Jacobian matrix. 



111 



with J{2D) — /i|i/2|2 — /i|2/2|i, we can express the single component of the curl and 
the divergence of g in terms of Lagrangian derivatives of g{X, t) and f{X,t): 

92,1 - 9i,2 = (5'2|i/2|2 - 5'2|2/2|i + fi'i|i/i|2 - 5'i|2/i|i) J'^ , (4a) 

91,1 + 92,2 = (5l|l/2|2 - 5l|2/2|l - 52|l/l|2 + 52|2/l|l) ■ (4b) 

Using the exact integrals (2d) and (2e) of eqs. (la) and (lb), eqs. (Ic) and (Id) are 
transformed into a closed Newtonian system for the trajectories in 2D: 

/2|l/2|2 — /2|2/2|1 + /l|l/l|2 — /l|2/l|l = , (5a) 
/l|l/2|2 - /l|2/2|l - /2|l/l|2 + /2|2/l|l = AJ{2D) - AttGq . (5b) 

(Notice that we have multiplied with J, i.e., we must ensure that J 7^ 0.) Looking at 
the equations (5) we see that we have just written out functional determinants: 

djfiJi) d{h,h) _ 

dUij2) Oihji) _^4(Z^^4.G^. (5b) 



d{Xi,X2) d{x,,X2) a(Xi,X2) 

In 3D we have to employ some more tensor algebra, but the procedure is the same. 
The relevant formula for the transformation of any tensor (here exemplified for the 
acceleration gradient {gi,j)) reads: 

9i,j 9i\khk,j — '^y^kim^jpq9i\kfp\lfq\m i (6) 



t We are talking about a purely two-dimensional space and not about cylinders in the third direction 
in which case we would have /s = X^. 
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where we have used the formula for the inverse Jacobian (2h) . In view of the definition 

— * — * — * 

of a functional determinant of any three functions A{X,t), B{X,t) and C{X,t), 
J{A,B,C) :— ' ' J = ek£mA kB.iC m , 

C(^Ai,A2,A3j 

e.g., for the Jacobian determinant we simply have J = <J{fi, f2, fs), we can write 
equation (6) as 

9i,j — '2j^3P<l'^{9ii fpi fq) ■ (6) 

The curl and the divergence of the acceleration field can be read from eq. (6) as the 
anti-symmetric part of the acceleration gradient and its trace (here, repeated indices 
imply summation as before, but with j, k running through the cyclic permutations of 
1,2,3): 

1 1 •■ 

9[i,3] = -2 ^ 9)k = i^^pq\jJUi]Jpi fq)J~^ , (7a,b,c) 

9i,i = (V • ^) = ^eabc J(L fb, fc)J-' . (7d) 

Inserting (7) and the exact integrals (2d) and (2c) into (Ic) and (Id) wc finally obtain 
the Lagrange-Newton system [7] (no background source, in particular A = 0) and 
[8] (including backgrounds of Friedmann type): 

J(/l,/l,/3)+J(/2,/2,/3) =0 , (8a) 

J(/2,/2,/l)+J(/3,/3,/l) =0 , (8b) 

J{flJlj2)+J{kf3j2) =0 , (8C) 
Jik /2, /3) + Jik /3, /l) + Jik fl, /2) 

-AWi,/2,/3) ^-^nGg . (8d) 

I have written this system here in explicit form to make working with these equations 
most convenient. (Alternative forms of these equations may be found in [21]; in that 
paper we also employ the calculus of differential forms, which makes the above derivation 
even simpler.) 

The Lagrange-Newton system (8) is the basic system of equations we want to study. 
Although these equations look, and in fact are more complicated than their Eulerian 
counterparts, they have proven to be as useful for finding exact solutions as well as 
perturbative approximations beyond the linear regime. Here, I note that the rules 
of determinant manipulations apply to these equations, and we may evaluate many 
problems analytically. 



Let us summarize the main conclusions of this section: 
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• The Eulerian evolution equations for the velocity and density fields (la) and (lb) 
are integrated exactly in the Lagrangian picture by (2d) and (2e), respectively. 

• The transformation of the Eulerian field equations (lc,d) yields a system of 
Lagrangian field equations for the acceleration field. 

• Using the integrals for the acceleration and density fields (2d) and (2e) 
we arrive at the Lagrange-Newton svstcm (8) which is a closed system of 
Lagrangian evolution equations for the field of trajectories . Density and velocity are 
no longer dynamical variables, but are replaced by the single dynamical field f{X,t). 

• The Lagrange-Newton system (8) is equivalent to the Euler-Newton system (1) as 

— * 

long as the mapping / is invertible ( J > 0) and, in particular, non-singular (J ^ 0); (for 
a proof see [21]). Note that the Lagrange-Newton system remains regular at caustics 
(J = 0) where the Eulerian representation breaks down. Whether the Lagrangian 
equations still describe the physics correctly in the regime J < will be discussed in 
Section 3. 

2. Lagrangian Dynamics of Eulerian Fields 

The question how to obtain the evolution of Eulerian fields from the Lagrangian 
description is easily answered: Knowing the Lagrangian solution of the transformation 
(2a), X — f{X,t), we have to express the Eulerian fields first in terms of /, and then 
use the inverse of this transformation (2g), X — h{x, t), to map the variable X back to 
Eulerian space, e.g., for the velocity we have: 

v{x,t) = 'f{h{x,t),t) . (9) 

The power of a Lagrangian description mainly relies on this implicit determination of 
the evolution of Eulerian fields. However, we can only come back to Eulerian space 

— * 

as long as h exists. In general, the inverse transformation can be multivalued (see the 
discussion in Section 3). 

The same rule applies to any Eulerian field, so we only have to find the corresponding 
formula which expresses it in terms of /. 

Let us now discuss another way of describing the evolution of Eulerian fields along 
trajectories. This will lead us to equations which are frequently discussed in the 
literature. These equations are useful to illuminate the power of the Lagrangian 
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formalism outlined in Section 1, but they do not determine the evolution of Eulerian 
fields per se, as will be shown below. 

We may ask for evolution equations which involve the Lagrangian time-derivative 
(2c) instead of the Eulerian time-derivative dt as in (la) and (lb). They may be written 
in symbolic form as: 

i.e., a system of u equations which determines the evolution of the u variables along 
the flow lines. A solution of this system of equations does not provide the full answer, 
which needs knowledge of the trajectories themselves! Wc will see that the procedure 
outlined at the beginning of this subsection does provide the full answer. We shall 
derive now Lagrangian evolution equations for a number of fields of interest, and shall 
then reinforce the Lagrange-Newton system (8). 

We start with eq. (la) and write it down in index notation, 

dt Vi + VkVi^k = 9i ■ (10a) 
Performing the spatial Eulerian derivative of this equation and using (2c) we obtain: 

= -%kVk,j + gi,j ■ (106) 

Eq. (10b) is an evolution equation for the velocity gradient (vij) along the flow hues 

— * 

/. It is convenient for the discussion of fluid motions to split it into its symmetric part 
(expansion tensor %), its antisymmetric part (vorticity tensor Wij), and to separate the 
symmetric part into a tracefree part (the shear tensor a^) and the trace (the rate of 
expansion) 9 Vi/. 
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u 



where v^ij) = ^{vij + vj^i) and = ^{vij - Vj,i). 

Inserting (11) into (10b) we obtain the following evolution equations (compare [20], [22], 
[32], [28] (§22): 

e = -^e^ + 2(Lo''-a^)+g,,i , (12a) 
2 

i^ij)' — —-^^^ij — CTikUJkj — ujikf^kj + 9[i,j] , (12b) 
2 2 1 

i^ij)' = ~'^^'^i3 ~ '^ik'^kj " ^^ik^kj '^('^ ~ ^ + 9(.i,j) — -^9k,k^ij , (12c) 

where :— \(Jij(Tij and o;^ := ^UijUiij. 

We can read these equations in the sense that they reconstruct the acceleration 
gradient {gi,j) in terms of kinematical variables: eq. (12a) gives the trace of {gi,j)-i eq. 
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(12b) its antisymmetric part, and eq. (12c) its tracefree symmetric part, the Newtonian 
tidal tensor E^j := g^^j - \gk,k5ij. 

We may now establish a system of equations for the variables 9, Uij and ujij by 
using the Euler-Newton system (1). This system only constrains the trace and the 
antisymmetric part of {gi,j)i but not the tidal tensor Eij] we arrive at: 

Q — —q9 , (13a) 

9 = -\9^ + 2(^2 _ ^2^) ^ ^ _ ^^Q^ ^ ^^3^^ 
2 

{ujiY ^ --9uji + (TijUjj , (13c) 

2 2 
{(Tij)' = --^Oaij - aikCFkj - uJikUJkj + ^(^ ~^ + ^ij ; (13d) 

(we have expressed the vorticity tensor in terms of the vector Cu — x v by means of 
the formula Uij = —eijkUJk)- 

There have been efforts in the literature to close this system of ordinary differential 
equations by using the corresponding equations of General Relativity (compare [1] and 
the lecture by Matarrese) and looking at their Newtonian limits. (In fact, the eqs. 
(13a-d) are formally identical to their GR counterparts in "comoving", i.e., Lagrangian 
coordinates.) However, it turns out that, although we can formally derive an evolution 
equation for the tidal tensor, the system of equations (13) supplemented by the evolution 
equation for the tidal tensor is not a system of ordinary differential equations and the 
problem remains "non-local" [2], [24], [23], [21]. 

As a matter of fact, in Newtonian theory we don't need evolution equations for 
tracefree symmetric tensors (like the shear and tidal tensors) to get a closed system 
of equations. We have already obtained such a system without them: the Lagrange- 
Newton system (8), which is a set of partial differential equations. 
To see the relation of the system of equations (13) to the Lagrange-Newton system we 
can show the following (the proof I leave to the reader as an excercise): 

2.) Eq. (13a), the continuity equation, is equivalent to eq. (lb) and is integrated in the 
Lagrangian framework by (2e). 

zz.) Eq. (13b), Raychaudhuri's equation, is equivalent to eq. (Id), if g is related to the 
velocity as in cq. (la). 

III.) Eq. (13c), the Kelvin-Helmholtz vorticity transport equation, is equivalent to eq. 
(Ic), if g is related to the velocity as in eq. (la). 

Note that eq. (13c) can also be integrated exactly in the Lagrangian picture, a result 
due to Cauchy (see [31] - a good textbook on Lagrangian dynamics): 

uj^ju;{X,to)-Vof . (14) 
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We arrive at the following conclusions of this section: 

• The equations (13a-c) are equivalent to the equations (Ib-d) provided the relation 
between velocity and acceleration is given by (la). 

• The equations (13b, c), if expressed in terms of Lagrangian coordinates, yield a 
closed Lagrangian svstem , the Lagrange-Newton system (8a-d). 

• The evolution of the tracefree symmetric tensors and E^j can be calculated 
after a solution to the Lagrange-Newton system is obtained, the formulas can be read 
off from eq. (6) (and a similar equation for the velocity gradient): 

(^ij ~ '^^^jpq>^ ifii fpi fq) gj'^"^'^'^ i.foi fpi fq)^ij 1 (15a) 
Eij = '^y^jpq'^ifii fpi fq) '^J^opq'^ifoi fpi fq)^ij (l^^) 
= ^^OPqJ' ifii fpi fq) - ^ - ■ ^^^^"^ 

We could use eq. (15c) to give another way of stating the Lagrange-Newton system: it 
is equivalent to the conditions that Eij is symmetric and tracefree (compare (7)): 

^[i J] = ^ (8a, 6, c) ; 
Eu^Q ^ {8d) . 



• The formulas (14) and (15) arc examples of formulas which express some Eulcrian 
field quantity in terms of /. Given any solution for /, we can insert it and use the 
inverse of the same solution to map the field into Eulerian space as in (9). (Integrals 
of the Eulerian evolution equations along the trajectories without knowledge of the 
trajectories themselves do not provide this information directly.) The perturbation 
solutions discussed in the next section will provide examples. 



3. Lagrangian Perturbation Theory 

The previous sections have equipped us with the necessary framework in which the 
Lagrangian dynamics of self-gravitating flows can be studied. We have reduced the 
description of the dynamics of any Eulerian field to the problem of finding the field of 

— * 

trajectories / as a solution of the Lagrange-Newton system (8). This problem will be 
addressed now. 
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3.1. The Lagrangian perturbation approach 

As in the Eulerian case we are not able to write down a solution for generic initial data 
(i.e. without any symmetry assumptions like plane or spherical symmetry). We may 
start with the simplest class of solutions, the homogeneous-isotropic ones, and then 
investigate a perturbative treatment of inhomogeneities. For homogeneous solutions 

— * 

the deformation tensor {fi^k) does not depend on X and the trajectory field is that of a 
homogeneous deformation: 

fi^{X,t) = aij{t)Xj ; aij{to) := 6ij ; (16a, 6) 

we assume hereafter that it is isotropic, 

aij{t) = a{t)6ij . (16c) 

Inserting this ansatz into the Lagrange-Newton system (8) yields the well-known 
equation 

a ^ A - A^Gen ^^^^^ 
a 3 

where qh — QnitojaT^ is the homogeneous density field calculated from (16a,c) and 
(2e). We may use it to integrate eq. (16d) yielding Priedmann's differential equation: 

a? -e SttGqh + A 

— = ; e = const . (loej 

a"' 3 

Solutions of (16e) cover the standard Friedmann cosmologies. We now define the 
(not necessarily small) deviation p from the "background model" f^, so that the full 

— * 

displacement map / of an inhomogeneous model reads: 

f^f"+p{X,t) ; p{X,to):^0 . (17a, 6) 

(Eq. (17b) expresses the freedom of starting at the same initial position as in the 
homogeneous case.) 

The Lagrangian perturbation approach then consists of solving the Lagrange-Newton 
system for powers of e to obtain the mth order solution: 

m 
i=l 

Note the following important differences to the Eulerian perturbation approach: 1. we 

— * 

need not perturb the density field as in the Eulerian case; 2. the perturbed fiow / 
is expressed in terms of coordinates which are constant along this flow, while in the 
Eulerian case the perturbation is expressed in a flxed frame, which does not change 
by moving away from the homogeneous flow. This implies that, although deviations p 




Figure 2. The unperturbed Hubble flow is represented by the motion of two fluid 
particles (dashed lines) and the perturbed flow is superimposed (full lines), Fig. 2a. In 
spite of the smallness of the deviations from the Hubble~flow compared to the Hubble- 
displacements, the density in Eulerian space can experience large changes, Fig. 2b. 



from / might be small compared to the homogeneous displacement, the Eulerian fields 
evaluated along the perturbed flow can experience large changes (Fig. 2). 

Solutions of the form (17) can be found and in the following subsections we shall 
learn more about them. 



3.2. A comment on average flows 

In the ansatz (17) we have splitted the inhomogeneous model into a homogeneous- 
isotropic model and a deviation from it. Before we discuss solutions of the form (17) 
we have to look at another assumption which has to be imposed on the inhomogeneous 
models. This will be necessary for their application to cosmology. 

In cosmology we want the homogeneous-isotropic deformation to describe the mean 
motion of the whole universe, and we want to interprete this motion as a reference 
frame (associated with, e.g., the microwave background). This interpretation requires 
that the background model (which, so far, is a homogeneous solution) correctly describes 
the average flow, i.e., inhomogeneities or peculiar-motions must vanish on average on 
some large scale. This we call the cosmological principle of homogeneity. That the 
average flow is isotropic is another assumption which we set a priori. It is, however, 
supported by the extreme isotropy of the microwave background. Therefore, in order 
to make Lagrangian perturbation models applicable to reality, we have to carefully 
analyze the average of inhomogeneous models. For this purpose we may look at the 
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general expansion law in Newtonian cosmology by averaging the system of equations 
(13). A general discussion is to be found in [15]. This study shows that the general 
expansion law for an expansion function ax> '■= V^^^, where V is the volume of the 
spatial domain ©(i) on which the average is performed, is of the form: 

3— +47rG'^-A = fap3 / V ■ u] +a:^^ [ V-Q , Q -.^ (uV ■ u - u -Vu) , 

av 3 \ Jv{t) ) Jv{t) 

(18a) 

where My is the total mass contained in V{t). Eq. (18a) shows that the average 
expansion is of Priedmann type, av = a, only in the case where the averages of the 
nonlinear terms on the r.h.s. (which depend on the inhomogeneities) vanish (compare 
(16d)). In other words, if we follow inhomogeneities into the nonlinear regime (and that's 
what we have in mind with Lagrangian approximations), then we actually neglect the 
influence of the inhomogeneities on the global expansion when we use homogeneous 
solutions for the average flow. Since these nonlinear source terms can be written as a 
divergence [15], we can use GauB's theorem to transform the volume integrals in (18a) 
into surface integrals over the boundary of the domain &D, e.g.: 

/ d^xV-Q^ [ dS-Q . (18b) 

Jv{t) JdV{t) 

— * 

The flux of u and Q through the boundary of the averaging volume is assumed to 
be negligible. This can only be true, if the extent of structures in the Universe is 
considerably smaller than the averaging domain. In standard cosmologies we assume 
that this is indeed the case on some large scale where we apply our models, and we even 
assume that this boundary flux is exactly zero. 

It is interesting to see at this stage that we implicitly satisfy this assumption by 
imposing periodic boundary conditions on the cosmic peculiar-fields: Periodicity means 
that we can alternatively view the Universe as a torus. A torus is a compact space 
without boundary, and consequently any flux through the boundary vanishes exactly. 

3.3. Perturbation solutions at first order 

In this subsection wc shall write down the general first-order solution in order to see 
that we need quite a list of arguments to reduce it to what we call the "Zel'dovich- 
approximation" [34], [35]. 

A solution to the basic system of equations (8) is general, if we can specify 
initial data for the velocity independently of a given density field (or acceleration 
field, respectively). If we split these data into their homogeneous parts (belonging 
to the background solution) and deviations from the background, peculiar-velocity u 
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and peculiar-acceleration w, then we have to specify 

U :^u{X,to) =:U^ + ; W:=w{X,to) , (19a) 

— * _ 

where we have sphtted the pecuhar-velocity into its irrotational (longitudinal) part U 
and its rotational (transverse) part C/^ (the acceleration is irrotational in Newton's 
theory) . We may introduce scalar potentials and a vector potential such that 

=: Vo-S ; =:VoxS ; W =: Vo$ • (196) 

Solving the Lagrange-Newton system to the first order with the perturbation ansatz 
(17) we obtain [9]: 

f'^ = a{t)X + b^{t)V,S^^\X) + 62(t)Vo X S'^^\X) + h{t)Vo^'^^\X) , (20) 

where the bi{t) depend on the chosen background solution a{t), bi{to) = 0, and the X- 
dependent perturbation potentials have to obey Poisson equations which relate them to 
the initial data (19b): 

Ao^^') = AoSto ; AoS^'^ = Aq^^o ; Aq^^'^ = Ao^l . (20a - e) 

(Aq denotes the Laplace operator with respect to Lagrangian coordinates.) 
The solution (20) shows that the evolution of a fluid element depends (non-locally) 
on the distribution of all fluid elements, because Poisson equations have to be solved 
(which is clear in view of the structure of the theory). The rotational part stems from 
the Lagrangian equations (8a-c) which must not be forgotten! Only at flrst order the 
rotational part directly represents rotational flow in Eulerian space. At higher orders 
rotational parts arise also for irrotational flow. This can be explained by the following 
consideration: Imagine an arrow for which we require that it keeps a fixed angle to some 
axis in an Eulerian coordinate system. The arrow may be placed on a system which is 
orbiting around some point. In this (Lagrangian) system the arrow performs a rotating 
movement, while in Eulerian space it is just translationally shifted. Indeed, if we restrict 
the motion to be irrotational in Lagrangian space, then this restricts possible classes of 
motion far more than the requirement of irrotational motion in Eulerian space. 

Let us return to the solution (20). At first order it is straightforward to solve the 5 
Poisson equations (20a-e) exphcitly: 

^(1) = + * ; s^^'' =sto + ^ ; = ml + n , 

Ao* = ; Ao^ = ; AqQ = . (21a - e) 

— * — * 

We can get rid of the harmonic functions and Q and we can make -S", 5" and $ 

unique, if we impose periodic boundary conditions on the initial data: the only periodic 
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harmonic functions are constant in space, a constant which expresses a translation of 
the whole simulation box and which can be set to zero, since the basic equations are 
translationally invariant (see [12], Appendix C and [21]). If we don't do that then the 
solutions are not unique! Removing this freedom reduces (20) to a local approximation 
which is only possible at first order. 

For the purpose of simulating the evolution of large-scale structure forward in time 
we can simplify the solution further by assuming that, initially 

U{X) = W{X)to , (22) 

i.e., peculiar-velocity and -acceleration are parallel to start with (in particular, we 
have U'^{X) — 0). This assumption is justified on physical grounds: as long as 
the perturbations are small, we may describe them by the (Eulerian) linear theory. 
Then, asymptotically (after the decaying modes have died out) the condition (22) is 
approximately reached and we may start with this condition from the outset. Employing 
(22) the first-order solution reduces to the "Zel'dovich-approximation" (compare [9] for 
the Einstein-de Sitter background, and [8], [4] for general backgrounds including a 
cosmological constant to see that the time-coefficients come out right): 

fl^^ ^a{t)X + b{t)VoS{X) . (23) 

Note that, strictly speaking, (23) is not the growing part of the general solution, but 
is that part which results from the restriction (22) on initial data: this restriction 
automatically cancels the decaying modes at first order, but not at higher orders. 



3.4- Higher-order solutions 

The solution of the Lagrange-Newton system for higher orders in e for the initial data 
(22) systematically gives us higher-order corrections to the "Zel'dovich-approximation" . 
For general initial data, already the second-order solution needs four pages to show it 
(compare [11] for the longitudinal part, which comprises the solution of one of the four 
equations (8a-d)). Restricting them to the setting (22) we obtain the much simpler 
form which is commonly being studied ([5], compare Bouchet's lecture, and [10], [11] 
for alternative models of this form) : 

/t^l = a{t)X + q,{t) Vo5W(X) + q,{t) VoS^'\X) . (24) 

The perturbation potentials in (24) have to be constructed by solving iteratively the 
two Poisson equations: 

Ao^W = /(^i,|fc)io , (24a) 
Ao5(^) = 277(5« ) . (24b) 



16 



/ and // denote the first and second principal scalar invariants of the tensor gradient 
in brackets: 

I{S^i^k)^tr{S^i^k)^^oS , (24c) 

iHSll^,) = llMSll^,))' - HiSll^.n . (24d) 

The third-order solution for the setting (22) may be found in [12], the time-coefficients 
being evaluated for an Einstein-dc Sitter background. (For other backgrounds we cannot 
write them in explicit form, see [18] for a thorough discussion.) Illustrations are given 
in Figs. 3,4. 



3. 5. Comparison of Eulerian and Lagrangian perturbation theory 

For the comparison of both perturbation approaches we look at the equations 
governing the evolution of inhomogeneitics in the first-order Eulerian and Lagrangian 
approximations. Consider the contrast density A := (f? — Qh)IQ) —oo < A < 1, 
which is more adapted to the nonlinear situation than the conventional definition 
S — {g — Qh)/qh = A/(l — A), defined in Eulerian perturbation theory. For this 
field we can find the following exact evolution equations by inserting the continuity 
equation into Raychaudhuri's equation (13b) and splitting off the homogeneous part 
[29], [8], [9]: 

A = (A - 1)/ , (25a) 

A + 2-A - 47rG'^jfA = (A - 1)277 , (25b) 

a 

where I and 77 denote here the invariants of the peculiar-velocity tensor gradient with 
respect to Eulerian coordinates which are scaled by the expansion factor, qi = Xi/a{t): 
(dui/dqj). 

For II — the equations (25) are (except for the Lagrangian time derivative and 
the redefinition of the density contrast) the same equations as known in Eulerian linear 
theory for S{q,t) [28]: 

I 5^ = -f , (25a)^ 

5' + 2— 5'-At:Gqh5' = Q . (25b)^ 
1 adt " 



d1 



1 



These equations can be obtained by linearizing (25). On the contrary, the Lagrangian 
hnear solution solves the exact equations (25) for 77 = 0. 

I emphasize that the similarity between the linear and the nonlinear case restricted 
to 77 = is nontrivial: a nice excercise is to compute the equation for the conventional 
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N= 1793, dt=0.40E-03, eps=0.50E-01, toll=0.50E+00 fc=2.0 



Figure 3. The trajectories for a simple plane wave model (2D) are plotted for the first- 
order (left, top panel) and for the second-order solution (left, bottom panel) (taken 
from [11]) compared with the trajectories in a two-dimensional tree-code simulation 
(taken from [17]). Those trajectories are plotted which end in an Eulorian spacctime 
section to show that the second-order model displays a secondary shell-crossing (this 
is not visible for trajectories which lie in this plane for all times - they follow a 
plane-symmetric collapse for which the first-order model is the exact solution). Two 
generations of caustics successively appear, whereas the numerical simulation shows 
further generations. 

density contrast S from the nonlinear equation (25b) for 11 = (solved by the 
Lagrangian first-order solution) using the definitions A := S/{l+S) and ' := dt\q+^u-'Vq, 
and to compare with the Eulerian linear equation (25b)^: 

S + 2-6 - AttCphS + 
a 



65 - 26^ + 2-SS - AttGphS'' = 
a 



(26) 
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Figure 4. A slice tlirougli the density fields (512'^ trajectories collected into a 256"^ 
pixel grid) predicted by the Lagrangian approximations at first order (left) and third 
order (right panel) for a 3D plane wave model is shown (taken from [17]). The second- 
and third-order schemes feature secondary shell-crossings. This hierarchy of nested 
caustics continues down to smaller scales as was shown in a two-dimensional numerical 
simulation in [19]. 



Hence, first-order Lagrangian perturbations involve nonlinearities due to nonlinearities 
in the dependent variable 5, but also due to products of u and 5 arising from the 
convective derivative (contained in the overdot). 

This exercise demonstrates the inherently nonlinear character of a Lagrangian 
perturbation approach. Essentially, this property can be traced back to the implicit 
solution of the Eulerian convection of the flow {u ■ Vq)u by the Lagrangian time- 
derivative, and to the fact that the density is exactly known for any solution. 

The second-order Lagrangian theory takes the term II in equation (25) 
approximately into account, and thus covers essential effects of the tidal action on the 
fluid. 



3.6. Limitations of the Lagrangian perturbation theory 

In order to recognize the limitations of the Lagrangian perturbation approach, two 
classes of arguments can be given: the flrst is analytical, based on the approach itself; 
the second is furnished by comparisons with N-body simulations. Since the latter is the 
topic of another lecture (by Peter Coles, and ref. therein), I here concentrate on the 
former, and only give two illustrations of the latter. 

Strictly speaking, the Lagrangian perturbation approach is limited by construction 
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to the regime t < tc„ where 

\Pi^j(Q\ «a(Q . (27) 

This condition is very conservative in view of the success of these approximations if 
followed up to shell-crossing, i.e., \pi,j{tc)\ — 0{a{tc)) and beyond (see Coles' lecture). 
The condition (27) has its counterpart in the condition 6 « 1 in Eulerian perturbation 
theory. At the epoch of shell-crossing, all orders of the perturbations yield displacements 
which are of comparable magnitude; consequently the perturbation approach breaks 
down there (we only expect it to converge before shell-crossing). In practice, however, 
we want to follow the solutions beyond shell-crossing (the Lagrange-Newton system is 
still regular in this regime, i.e., the individual trajectories don't feel their intersection 
with others) . This demands that we must argue whether and how this extrapolation is 
justified. The examples given in the previous subsection suggest that the extrapolation 
can indeed be meaningful. 

— * 

As mentioned earlier, in the regime after shell-crossing J < 0, the mapping / is no 

— * 

longer unique and the inverse h is multivalued. If we assume that the density can be 
written as a superposition of the partial densities corresponding to the different streams, 
then a reasonable extension into the multi-stream regime for the density is 

g[x, t] = E ^(^7(^' t),t) = J2 ^'^^'"'i'"^^^ ; (28) 

7=1 7=1 |det(/j|fe(/i^(f,t),t) )| 

the total density in the n-stream region is the sum of the moduli of the individual 
densities of the streams. This extension is imphcitly used when we construct the 
density field from displaced particles: the formula (28) is approximated by the collection 
of particles in an Eulerian grid at times after shell-crossing. However, due to the 
presence of self-gravity the streams interact and the formula (28) has to be taken with 
caution. Indeed, expressing the superposition assumption in terms of the accelerations 
we encounter a problem: In the multi-stream region there are n fluid elements at the 
same Eulerian position, thus the total gravitational field strength at this position would 
have to be added up similar to (28) and the fluid particle should be pushed by this 
total strength. On the contrary, the Lagrangian solutions still push the fluid elements 
with their individual accelerations. In other words, Einstein's equivalence principle of 
inertial and gravitational mass is violated: the particles are not accelerated by the total 
gravitational field strength. We have to find extensions of the theory presented (e.g., in 
the framework of the Vlasov-Poisson system) to properly deal with this problem. 

The performance of the Lagrangian perturbation solutions after shell-crossing as 
inferred from a number of comparisons with N-body simulations appears to be weak 
due to the shortcomings discussed above. This is most easily visible in pancake models 
[13]. If we consider hierarchical cosmogonies, especially models which involve much 
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small-scale power, then the power spectrum of fluctuations must be sufficiently steep to 
avoid shell-crossing on small scales. Otherwise, these high-frequency components have 
to be truncated or filtered (see Coles lecture and ref. therein). In the former 
spectral index n = —3 in a power law spectrum oc Ikl"" lies close to the limit between 
"pancake-type" and "hierarchical" cosmogonies having n > —3. This is so, because 
in the case n = —3 the integrated power per wave number interval is approximately 
constant (logarithmic), so that all waves enter the nonlinear regime at about the same 
time (Fig. 5). For more small-scale power the smoothing scale of high-frequency modes 
has to be close to (but smaller than) the nonlinearity scale. In Fig. 6, 1 present the result 
of a comparison for a variant of the CDM cosmogony, the BSI model (see Gottlober's 
lecture, this volume) with less small-scale power than standard CDM. 




Figure 5. A slice of tlie density field is shown for a PM simulation (left) and the 
second-order Lagrangian approximation (right panel) (taken from [27]). The spectrum 
of fluctuations was given by a powerlaw with index n = —3, which can be taken to 
discriminate between the "pancake" and "hierarchical" regimes of structure formation. 



Let US again summarize the main conclusions of this section: 

• In a Lagrangian perturbation approach only one variable, the trajectory field /, is 
perturbed. The perturbation is evaluated along the perturbed orbits. Although the 
perturbations of the orbits might be small, the Eulerian fields evaluated along the 
perturbed orbits can experience large changes. 

• The first-order Lagrangian approximation reduces to the "Zel'dovich-approximation" . 
if the initial data are restricted such that peculiar-velocity and peculiar-acceleration are 
parallel. 
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Figure 6. A slice of the density field in a 200h~-'^Mpc box is shown for a PM 
simulation (left) and the second order Lagrangian approximation (right panel) (taken 
from [33]). The spectrum of fluctuations was a (COBE normalized) non-standard 
CDM spectrum with broken scale invariance (BSI). The spectrum had to be smoothed 
at the high-frequency end for the Lagrangian scheme to avoid substantial post- 
singularity evolution. 

• Imposing periodic boundary conditions on the cosmic peculiar-fields is necessary to 
guarantee 1. uniqueness of the solutions, and 2., to assure that the average flow is given 
by the homogeneous "background" solutions. 

• Higher-order effects amount to additional internal structures such as "second generation" 
pancakes, filaments and clusters. They also accelerate the collapse. 

• A Lagrangian perturbation solution inherently includes nonlinear terms. This is also 
true for the first-order solution, which is linear in Lagrangian space. 

• The temporal limit of application of the Lagrangian schemes is reached at the epoch 
of shell-crossing, where the displacements for each order attain similar magnitudes. 
Following these schemes beyond shell-crossing time requires truncation of high- 
frequency modes in the initial fiuctuation spectrum, thus defining a spatial lower limit 
of application, which roughly corresponds to galaxy group mass scales, lO^^M© [26], 
[16]. 

The Lagrangian perturbation approach is investigated to arbitrary order in [21]. In that 
paper as well as in the review articles [3], [6], [14] and [30] you may find references to 
related work. 
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